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ABSTRACT 

Numerical simulations indicate that the smooth, photoionized intergalactic medium 
(IGM) responsible for the low column density Lya forest follows a well defined 
temperature-density relation, which is well described by a power-law T = To(p/p)' y ~ 1 . 
We demonstrate that such an equation of state results in a power-law cutoff in the 
distribution of line widths (6-parameters) as a function of column density (iV) for the 
low column density (N ^ 10 145 cm -2 ) absorption lines. This explains the existence 
of the lower envelope which is clearly seen in scatter plots of the 6(iV)-distribution 
in observed QSO spectra. Even a strict power-law equation of state will not result in 
an absolute cutoff because of line blending and contamination by unidentified metal 
lines. We develop an algorithm to determine the cutoff, which is insensitive to these 
narrow lines. We show that the parameters of the cutoff in the 6(iV)-distribution are 
strongly correlated with the parameters of the underlying equation of state. We use 
simulations to determine these relations, which can then be applied to the observed 
cutoff in the &(iV)-distribution to measure the equation of state of the IGM. We show 
that systematics which change the 6(7V)-distribution, such as cosmology (for a fixed 
equation of state), peculiar velocities, the intensity of the ionizing background radia- 
tion and variations in the signal to noise ratio do not affect the measured cutoff. We 
argue that physical processes that have not been incorporated in the simulations, e.g. 
feedback from star formation, are unlikely to affect the results. Using Monte Carlo 
simulations of Keck spectra at z = 3, we show that determining the slope of the equa- 
tion of state will be difficult, but that the amplitude can be determined to within ten 
per cent, even from a single QSO spectrum. Measuring the evolution of the equation 
of state with redshift will allow us to put tight constraints on the rcionization history 
of the universe. 

Key words: cosmology: theory - intergalactic medium - equation of state - hydro- 
dynamics - quasars: absorption lines 



1 INTRODUCTION 

The intergalactic medium (IGM) at high redshift (z ~ 2-5) 
manifests itself observationally by absorbing light from dis- 
tant quasars. Resonant Lya absorption by neutral hydrogen 
along the line of sight to a quasar results in a fluctuating 



Lya trs nsmission (optical depth) . Regions of enhanced den- 



evolution of the quasar population and the star formation 
history of the universe. 

Hydrodynamical simulations of structure formation in 
a universe dominated by cold dark matter and including 
an ionizing background, have been very successful in ex- 



plain i ng the properties of the Lya forest ( s.g. Cen et al 
f994|; Zhang, Anninos fc Norman 1995; Hcrnquist ct al. 199(: 



sity give rise to increased absorption and appear as a forest 



Miralda-Escude et al. 1996; Theuns, Leonard & Efstathiou 



199^ ; [Theuns et al. 1996] ; pave et al. I999| ). They show that 
the low column density (N % 10 1 *' cm"') absorption lines 
arise in a smoothly varying low-density (8 ^ 10) IGM. Since 
the overdensity is only mildly non-linear, the physical pro- 
cesses governing this medium are well understood and rela- 
tively easy to model. On large scales the dynamics are de- 



of abso 'ption lines bluewards of the quasar's Lya emission 



line. The fact that not all the light is absorbed imp lies that 
the IGM is highly ionized ( Gunn fe Peterson 1965 ), but the 



time and origin of the reionization of the gas are still un- 
known. Since quasars and young stars are sources of ionizing 
radiation, the ionization history of the gas depends on the 
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termined by gravity, while on small scales gas pressure is 
important. The existence of a simple physical framework 
and the abundance of superb observational data make the 
Lya forest a n ext remely promising cosmological laboratory 
(see Rauch (1998) for a review). In this paper, we shall in- 
vestigate the effect of the thermal state of the IGM on the 
forest. 

For the low-density gas responsible for the Lya forest, 
shock heating is not important and the gas follows a well 
defined temperature-density relation. The competition be- 
tween nhotoionization heating and adiabatic cooling results 



0.025) significantly broadens the absorption lines, although 
some discrepancy with observations may remain. One way 
to increase the temperature further would be to change 
the reionization history. Uncertainties in the redshift of He 
reionization, in particular, can affec t the temperature of the 
IGM. Haehnelt & Steinmetz (1998) demonstrated that dif- 



ferent reionization histories result in observable differences 
in the ^-distribution. Other mechanisms that have been pro- 
posed to boost the temperature are photoelect ric heating 
of dust grains ( Nath, Sethi fe Shchekinov 199S ), Compton 



in a po Ler-la w 'equation of state' T = Tn(p/ p) 1 (Hhii & 



heati ng by the hard X-ray background ( Madau fc Efstathiou 



199E ) and radiative transfer effects associated with the ion- 



jqq TJ) This equation of state depends on co smology ization of Hen by QSOs in the optically thick limit jAbel & 



iniy.a.tinn history For models with abrupt 



rniomza- 



Gnedin 
and rei 

tion, the IGM becomes nearly isothermal (7 fa 1) at the 
redshift of reionization. After reionization, the temperature 
at the mean density (To) decreases while the slope (7 — 1) 
increases because higher density regions undergo less expan- 
sion and increased photoheating. Eventually, when photo- 
heating balances adiabatic cooling as the universe expands, 
the imprints of the reionization history are washed out and 
the equation of state approaches an asymptotic state, 7 = 

1.62, To oc [fib/iVVOm/i 2 ] 171 ' 7 flHui fc Gnedin 1997p . Since 
the reionization history of the universe is still unknown, the 
physically reasonable ranges for the parameters of the equa- 
tion of state ar e very large (10 3 ° K < To < 10 4 5 K and 
1.2 < 7 < 1.7 flHui, Gnedin fe Zhang 1997|; [Hui fc Gnedin 



Haehnelt 1999) 



1997)) 



The smoothly varying IGM gives rise to a fluctuating 
optical depth in redshift space. Many of the optical depth 
maxima can be fitted quite accurately with Voigt profiles. 
The distribution of line widths depends on the initial power 
spectrum, the peculiar velocity gradients around the density 
peaks and on the temperature of the IGM. However, there 
is a lower limit to how narrow the absorption lines can be. 
Indeed, the optical depth w ill be smoothed on a sc ale deter- 
mined by three processes ( |Hui fc Rutledge 1997 ): thermal 
broadening, baryon (Jeans) smoothing and possibly instru- 
mental, or in the case of simulations, numerical resolution. 



Unfortunately, the fa-distribution is not very well suited 
for investigating the thermal state of the IGM. Although 
some of the broad lines correspond to density fluctuations 
on scales th at are affected, amo ng other things, by thermal 
smoothing (Theuns et al. 1999b), many are caused by heavy 
line blending and continuum fitting errors. The cutoff in the 
fa-distribution, on the other hand, will depend mainly on the 
temperature of the IGM and is therefore potentially a pow- 
erful statistic. In practice its usefulness is limited because 
many narrow lines occur in the wings of broader lines. Such 
narrow lines are often introduced by numerical Voigt profi le 
fitting algorithms (such as VPFIT flCarswell et al. 1987 )) 
to improve the quality of the overall fit to the quasar spec- 
trum. If the physical structure responsible for the absorption 
does not consist of discrete clouds, then the widths of these 
blended lines will have no relation to the thermal state of the 
gas. Furthermore, the number and widths of these narrow 
blended lines depend on the Voigt profile fitting algorithm 
that is used and on the signal to noise of the quasar spec- 
trum. 

The fa-parameter distribution is usually integrated over 
a certain column density range. Since the fa-distribution 
might depend on column density (N), more information is 
contained in the full fa(JV)-distribution. Scatter plots of the 
6(iV)-distribut ion have been pub li shed for many observed 



The firiiL Lwu depend un Lire Lhemial slale uf Lhe gas. While 
for hig h - rc3olution observations (cchcllc spectroscopy) the 
effective smoothing scale is not determined by the instru- 
mental resolution, numerical resolution has in fact been the 



QSO spectra (|e.g. Hu et al. 1995; 



Lu et al. 1991: ; Kirkman 



Tytler 1997} [Kim et al. 1997| ). These plots show a clear 



limiting factor in many simulations (see Theuns et al. ( 1998 ) 
for a discussion). 

The distribution of line widths is generally expressed as 
the distribution of the widths of Voigt profile fits to the ab- 
sorption lines, the fa-parameters. While the first numerical 
simulations showed good agreement with the observed fa- 
parameter distribution, higher resolution simulations of the 
standard cold dark matter mo del produced a larger fraction 
of narrow lines than obs erved (Theuns et al. 1998; Bryan et 
al. 1998j) . Theuns et al. ( 19981) suggested that an increase in 
the temperature of the IGM might broaden the absorption 
lines, while Bryan et al. ( 1998 ) argued that the most natural 



cutoff at low fa-parameters. However, this cutoff is not abso- 
lute. There are som e nar row lines, especially at lo w col umn 
densities. Lu et al. ( 1996 ) and Kirkman & Tytler ( 1997 ) use 
Monte Carlo simulations to show that many of these lines 
are caused by line blending and noise in the data. Some con- 
tamination from unidentified metal lines is also expected. 

The cutoff in the b(iV)-di stribu tion increases slightly 
with column density. Hu et al. (1995) conclude from Monte 



way to broaden the lines is to change the density distribu- 
tion directly. Note that increasing the temperature will also 
change the density distribution of the gas through increased 
baryon smoothing (Theuns, Schaye & Haehnelt 1999). 



Carlo simulations that this correlation is primarily an ar- 
tifact of the much larger number of lines at lower column 
density and the incr eased scatter in the fa-determinations. 
Kirkman & Tytler (1997) however, conclude from similar 
simulations that the correlation between the lowest fa- values 
and column density is a real physical effect. (Note that this 
corr elatio n is different from the one reported by Pettini et 
al. (1990). They found a general correlation between the 



Theuns et al. (1999a) showed that changing the cos- 



mology (lowering fi m from 1.0 to 0.3 and doubling Qi,h to 



6-parameter and column density of all lines. This correla- 
tion was later shown to be an artifact of the line selection 
and fitting procedure (Rauch et al. 199E ).) A lower envelope 
which increases with c olumn density has also been seen in 
numerical simulations (Zhang et al. 1997). 
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Table 1. 


Models simulated 






Model 


a L m 


n A 


tt b h 2 


h 


r H i 


S 


1 





0.0125 


0.5 


HM/2 


Sb 


1 





0.025 


0.5 


HM 





0.3 





0.0125 


0.65 


HM/2 


Ob 


0.3 





0.025 


0.65 


HM 


L 


0.3 


0.7 


0.0125 


0.65 


HM/2 


Lb 


0.3 


0.7 


0.025 


0.65 


HM 



heated by an imposed uniform background of UV-photons 
that ori ginate s from quasars, as computed by Haardt & 
Madau (1996). This flux is redshift dependent, due to the 



In this paper we shall demonstrate that the cutoff in 
the 6(iV)-distribution is determined by the equation of state 
of the low-density gas. Furthermore, we shall show that the 
cutoff can be determined robustly and is unaffected by sys- 
tematics like changes in cosmology (for a fixed equation of 
state) and can therefore be used to measure the equation of 
state of the IGM. 

We test our methods using smoothed-particle hydrody- 
namic (SPH) si mula t ions of the Lya forest as described by 
Theuns et al. (19981; 1999a). The parameters of the simu- 
lations are summarised in section H. The relation between 
the 6(AQ-cutoff and the equation of state is investigated in 
section fcA which forms the heart of the paper. Section H con- 
tains a detailed description of the procedure used to fit the 
cutoff in simulated Keck spectra. Systematic effects are dis- 
cussed in section g. In section g we test the procedure using 
Monte Carlo simulations. Finally, we summarise and discuss 
the main results in section bj. 



2 SIMULATIONS 

We have simulated six different cosmological models, char- 
acterised by their total matter density Q m , the value of 
the cosmological constant S7a, the rms of mass fluctua- 
tions in spheres of radius 8/t -1 Mpc, as, the baryon den- 
sity Slj,/i 2 and the present day value of the Hubble constant, 
Ho = lOO/i kms -1 Mpc -1 . The parameters of these models 
are summarised in Table hi. In addition to these models, we 
simulated a model that has the same parameters as model 
Ob, but with the He I and Hell heating rates artificially dou- 
bled. This may provide a qualitative model of the heating 
due to radiative transfer effec ts during the reionization of 
helium flAbel fc Haehnelt 1999"| ). We will call this model Ob- 
hot. The amplitude of the initial power spectrum is nor- 
malised to the observed abundance of galaxy clust ers a t 
z = 0, using the fits computed by Eke, Cole & Frenk (1996). 
We model the evolution of a periodic, cubic region of the 
universe of comoving size 2.5h~ 1 Mpc. 

The code us ed is adapted from the HYDRA code of 
Couchman et al. (|l995), which uses smooth particle hydro - 
dynamics (SPH) QLucy 1970|; |Gingold fc Monaghan 1977| ); 
see Theuns et al. (1999a; 1998J ) for details. These simula- 
tions use 64 3 particles of each species, so the SPH particle 
masses are 1.65 x 10 6 (fif,/i 2 /0-0125)(/i/0.5)" 3 M Q and the 
CDM particles are more massive by a factor f2cpM/^6- This 



resolution is sufficient to simulate line widths reliably (The- 



uns et al. 1998; Bryan et al. 199£, note that in the hotter 



evolution of the quasar luminosity function. The amplitude 
of the flux is indicated as 'HM' in the Tjjj column of Ta- 
ble [l] (where Fjj : is the H I ionization rate due to the ionizing 
background). For the low Q,bh 2 models, we have divided the 
ionizing flux by two, indicated as 'HM/2'. We do not impose 
thermal equilibrium but solve the rate equations to track the 
abundances of H I, H II and He I, He II and He ill. We assume 
a h elium abundance of Y — 0.24 by mass. See Theuns et 
al. jl998| ) for further details. 

At several output times we compute simulated spec- 
tra along 1200 random lines of sight through the simulation 
box. Each spectrum is convolved with a Gaussian with full 
width at half maximum of FWHM = 8 kms -1 , then resam- 
pled onto pixels of width 3 kms -1 to mimic the instrumental 
profile and characteristics of the HIRES spectrograph on the 
Keck telescope. We rescale the background flux in the anal- 
ysis stage such that the mean effective optical depth at a 
given redshift in all models is the same as for the Ob model. 
This model h as a mean absorpti on in good agreement with 
observations (Rauch et al. 1997) (Ob has f e a = 0.93, 0.33 
and 0.14 at z = 4, 3 and 2). Finally, we add to the flux in 
every pixel a Gaussian random signal with zero mean and 
standard deviation a = 0.02 to mimic noise. The spectra 
cover a small enough velocity range to b e fit by a flat contin - 
uum, as chosen by a simple procedure ( Theuns et al. 199£ ) 
described as follows. A low average continuum is assumed 
initially, then all pixels below and not within 1 a of this level 
are rejected and a new average flux level for the remaining 
pixels is computed. The last two steps are repeated until the 
average flux varies by less than 1%. This final average flux 
level is adopted as the fitted continuum and the spectrum is 
renormalised accordingly. The absorption features in these 
mock observations are then fitted with Voigt profiles using 



an automated version of VPFIT ( Jarswell et al. 1987 



simulations numerical convergence will be even better than 
in the cooler model S, which was investigated in detail by 
Theuns et al. ( |l99S| )). 

We assume that the IGM is ionized and photo- 



Although the simulated models have different equations 
of state, they cover only a limited part of the possible param- 
eter space (To, 7). The models were originally intended to in- 
vestigate the d ependence of QSO ab sorption line statistics 
on cosmology (Theuns et al. 1999a). Changing the reion- 
ization history can lead to very different values of To and 
7. To quantify the relation between the cutoff in the b(N)- 
distribution and the equation of state, it is necessary to in- 
clude models covering a wide range of To and 7. We therefore 
created models with particular values of To and 7 by impos- 
ing an equation of state on model Ob. This was done by 
moving the SPH particles in the temperature-density plane. 
The new models have the same three components (low- and 
high-density power-law equations of state and shocked gas) 
as the original model, but a different equation of state for the 
low-density gas. In particular, the intrinsic scatter around 
the power- law is left unchanged. 



3 THE 6(iV)-CUTOFF AND THE EQUATION 
OF STATE 

Fig. hi shows a contour plot of the mass-weighted distri- 
bution of fluid elements (SPH particles) from a numerical 
simulation in the temperature-density diagram. Noting that 
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-2-101234 
log (p b /<p b >) 

Figure 1. The mass- weighted distribution of fluid elements (SPH 
particles) in the temperature-density plane for model Ob-hot at 
z = 3. Three contour levels are plotted. The number density 
of fluid elements increases by an order of magnitude with each 
contour level. The dashed line is a power-law (least absolute de- 
viation) fit to the temperature-density relation of the gas with 
density in the range —0.7 < log(pt,/p(,) < 0.3. Clearly, the vast 
majority of the low-density gas follows this tight temperature- 
density relation. The dotted line shows the result of a similar fit 
for model S. 



the number density of fluid elements increases by an order 
of magnitude with each contour level, it is clear that the 
vast majority of the low-density gas (pb/pb Ss 10) follows 
a power- law equation of state (dashed line) . The two other 
components visible in Fig. hi are hot, shocked gas which can- 
not cool within a Hubble time and colder, high-density gas 
for which Hi and Hell line cooling is effective. 

In Fig. Ha we plot the 6(iV)-distribution for 800 random 
absorption lines taken from the spectra of model Ob-hot at 
redshift z — 3. A cutoff at low 6-values, which increases with 
column density, can clearly be seen. In Fig. tjb only those 
lines for which VPFIT gives formal errors in both b and N 
that are smaller than 50 per cent are plotted. This excludes 
most of the lines that have column density N ^ 10 12 ' 5 cm~ 2 
as well as many of the very narrow lines below the cut- 
off. Although the formal errors of the Voigt profile fit given 
by VPFIT have only limited physical significance, lines in 
blends tend to have large errors. Since the fe-parameters for 
blended lines can have values smaller than the minimum set 
by the thermal smoothing scale (i.e. thermal broadening and 
baryon smoothing), these lines will tend to smooth out any 
intrinsic cutoff. Removing the lines with the largest relative 
errors therefore results in a sharper cutoff. A smaller max- 
imum allowed error would result in the removal of many of 
the regular, isolated lines. 

The large number of data points plotted in Fig. Ha ex- 
cludes the possibility that the slope in the cutoff is due to the 
large decrease in the number of lines with column density or 
the in crease in scatter with decreasing column density (Hu 
et al. ( 199a ) reached this conclusion from analysing spectra 
that had about 250 absorption lines each). 

The 6(JV)-distribution for the colder model S is plotted 
in Fig. Pp. Clearly, the distribution cuts off at lower b- values. 
Let us assume that the absence of lines with low 6- values is 
due to the fact that there is a minimum line width set by 



the thermal state of the gas through the thermal broadening 
and/or baryon smoothing scales. Since the temperature of 
the low-density gas responsible for the Lya forest increases 
with density (Fig. hi) , we expect the minimum b- value to in- 
crease with column density, provided that the column den- 
sity correlates with the density of the absorber. 

To see whether this picture is correct, we need to in- 
vestigate the relation between the Voigt profile parameters 
N and b, and the density and temperature of the absorbing 
gas respectively. Peculiar velocities and thermal broadening 
make it difficult to identify the gas contributing to the op- 
tical depth at a particular point in redshift space. Further- 
more, the centre of the Voigt profile fit to an absorption line 
will often be offset from the point where the optical depth is 
maximum. We therefore need to define a temperature and 
a density which are smooth and take redshift space distor- 
tions into account. We choose to use optical depth weighted 
quantities: the density of a pixel in velocity space is the sum, 
weighted by optical depth, of the density at all the pixels in 
real space that contribute to the optical depth of that pixel 
in velocity space. We then define the density corresponding 
to an absorption line to be the optical depth weighted den- 
sity at the line centre. The temperature corresponding to 
absorption lines is defined similarly. 

In Fig. a the optical depth weighted density and tem- 
perature are plotted for two random lines of sight (dashed 
lines in the middle two panels), as well as the flux (without 
noise), and real space density, temperature and peculiar ve- 
locity. The dashed curves in the top panels are the Voigt pro- 
files fitted by VPFIT, vertical lines indicate the line centres. 
Absorption lines correspond to peaks in the (optically depth 
weighted) density and temperature, which are strongly cor- 
related. Although the blends indicated by arrows can be 
traced back to substructure in the peaks, their profiles are 
mainly determined by the density and temperature of the 
gas in the main peaks. 

In Fig. H the optical depth weighted gas density is plot- 
ted as a function of column density for the absorption lines 
plotted in Fig. Eh. There exists a tight correlation between 
these two quantities. Note that lines with column densities 
Ja 10 13 cm -2 correspond to local maxima in underdense re- 
gions. The optical depth weighted temperature is plotted 
against the 6-parameter in Fig. Ba. The result is a scat- 
ter plot with no apparent correlation. This is not surpris- 
ing since many absorbers will be intrinsically broader than 
the local thermal broadening scale. In order to test whether 
the cutoff in the &(A r )-distribution is a consequence of the 
existence of a minimum line width set by the thermal state 
of the gas, we need to look for a correlation between the 
temperature and 6-parameters of the lines near the cutoff. 

Determining a cutoff in an objective manner is nontriv- 
ial because of the existence of unphysically narrow lines in 
blends. We developed a fitting algorithm that is insensitive 
to these lines. This algorithm is described in the next sec- 
tion. Fig. H is a scatter plot of the lines with column density 
in the range 10 125 cm -2 < N < 10 145 cm -2 . The cutoff fit- 
ted to this distribution is also shown (solid line). The lines 
that are used in the final iteration of the fitting algorithm, 
i.e. lines that are close to the solid line in Fig. |fj| do in- 
deed display a tight correlation between the temperature 
and 6-parameter (Fig. Hb). The dashed line in Fig. Hb cor- 
responds to the thermal width, b — (2fcflT/m p ) 1 ' 2 , where 
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log N HI (cm" 2 ) 



log N H1 (cm" 2 ) 



log N HI (cm" 2 ) 



Figure 2. The fe(A r )-distribution for 800 random lines from model Ob-hot (a,b) and S (c), at z = 3. The position of each line is indicated 
by a cross. Errors are not displayed. In panels (b) and (c) only those lines are plotted for which VPFIT gives formal errors Ab/b < 0.5, 



AN/N < 0.5. The dashed line is the cutoff for the lines plotted in panel (b) over the range 10 125 cm 2 < N < 10 145 c 



Removing 



the lines with the largest relative errors in the Voigt profile parameters sharpens the cutoff. The b(7V)-distribution of the colder model S 
clearly cuts off at lower b-parameters than the 6(A r )-distribution of model Ob-hot. 




100 150 200 
v hub (km/s) 



300 



50 



100 150 200 
v hub (knVs) 



300 



Figure 3. Two random lines of sight through the box of model Ob-hot at z = 3. From top to bottom the solid curves are the flux, gas 
density, temperature and peculiar velocity respectively, all as a function of Hubble velocity. The dashed curves in the top panel are the 
Voigt profiles fitted by VPFIT, vertical lines indicate the line centres. The arrows point at the line centres of lines with a b-parameter 
smaller than the cutoff value of the Hui-Rutlcdge formula (equation 101). These lines have column densities N = 10 12,7 cm -2 (left panel) 
and N = 10 129 cm - 2 (right panel), and 6-parameters 13.1 kms -1 (left panel) and 13.6 kms -1 (right panel). They appear as diamonds 
in Fig. H. The dashed curves in the second and third panels are the optical depth weighted density and temperature respectively. 
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^fRp + >+ H * 
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0.8 1.0 1.2 1.4 1.6 
log b (km/s) 



1.0 1.2 1.4 1.6 
log b (km/s) 



2.0 



Figure 5. The optical depth weighted temperature at the line centres as a function of 6-parameter. The left panel contains all lines 
plotted in Fig. 0b. In the right panel only those lines are plotted which have fe-parameters within one mean absolute deviation of the 
power-law fit to the b(A r )-cutoff plotted as the solid line in Fig. rj. Only for the lines near the cutoff are the 6-parameters correlated with 
the temperature of the absorbing gas. The dashed line corresponds to the thermal width, b = (2fcsT/m p ) 1 ' . 




13 14 

log N HI (cm" 2 ) 



15 



Figure 4. The optical depth weighted gas density at the line 
centres as a function of neutral hydrogen column density. Data 
points correspond to the lines plotted in Fig. Ep. The column 
density of a Voigt profile fit to an absorption line is strongly 
correlated with the density of the absorbing gas. The outlier at 
log(p/p"f,) = 2.50 corresponds to a line in the wing of a highly 
saturated line. 



1 

^ 




12.5 



13.0 



13.5 



log N HI (cm 2 ) 



14.0 



14.5 



Figure 6. The 6(iV)-distribution for those lines from Fig. Eh 
that have column densities in the range 10 12 ' 5 cm~ 2 < N < 
10 14 ' 5 cm -2 . The solid line is the fitted cutoff for this set of lines, 
the dot-dashed line is the cutoff for the colder model S. The hor- 
izontal dashed line indicates the minimum 6-value of lines used 
in the fit, lines with smaller b-parameters are indicated by dia- 
monds. This minimum 6-value corresponds to the cutoff of the 
best-fitting Hui-Rutlcdge function (see section M). 



m p is the mass of a proton and fcs is the Boltzmann con- 
stant. Lines corresponding to density peaks whose width in 
velocity space is much smaller than the thermal width, have 
Voigt profiles with this fa-parameter. Since the temperature 
plotted in Fig. M is the smooth, optical depth weighted tem- 
perature, we do not expect the relation between T and fa to 
be identical to the one indicated by the dashed line, even 
if all the line widths were purely thermal. Although other 
definitions of the density and temperature are possible and 
will give slightly different results, qualitatively the results 
will be the same for any sensible definition of these physical 
quantities. Figs. W and 0b therefore suggest that the cutoff 



in the fe(7V)-distribution should be strongly correlated with 
the equation of state of the absorbing gas. 

Let us look in more detail at the relation between the 
fo(iV)-cutoff and the equation of state. We have already 
shown (Fig. |l|) that almost all the low-density gas follows 
a power-law equation of state: 

log(T) = log(To) + (7 - 1) log(p/p). (1) 

The relations between the density /temperature and column 
density/fa-parameter of the absorption lines near the cutoff 
can also be fitted by power-laws (Fig. ^ and Fig. |5p): 

log(p/p)=A + Blog(N/No), (2) 
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log(T) = C + Dlog(6). 



(3) 



Combining these equations, we find that the b(JV)-cutoff is 
also a power-law, 



log(b) = log(feo) + (T - 1) log(N/N ), 
whose coefficients are given by, 

log(bo) = i [log(To) - C + (7 - 1)4 , 



1-|(7- 



1). 



(4) 



(5) 

(6) 



Hence the intercept of the cutoff, logfeo, depends on both the 
amplitude and the slope of the equation of state, while the 
slope of the cutoff, V — 1, is proportional to the slope of the 
equation of state. If we set No equal to the column density 
corresponding to the mean gas density, then the coefficient 
A vanishes and log bo no longer depends on 7. 

The cutoff in the 6(A r )-distribution is measured over 
a certain column density range. We choose to measure the 
cutoff over the interval 10 125 cm" 2 < N < 10 145 cm" 2 . This 
interval corresponds roughly to the gas density range for 
which the equation of state is well fitted by a power-law. 
For lines with column density N ^ 10 12 ' 5 cm" 2 , the scatter 
in the observations becomes very large due to noise and line 
blending. The relation between the gas overdensity and H I 
column density depends on redshift. Hence different column 
density intervals should be used for different redshifts if one 
wants to compare the equation of state for the same density 
range. 

Fitting a cutoff to a finite number of lines introduces 
statistical uncertainty in the measured coefficients. We min- 
imize the correlation between the errors in the coefficients 
by subtracting the average abscissa value (i.e. log N) be- 
fore fitting the cutoff. This can be done by setting log No in 
equation |4j equal to the mean log N of the lines in the col- 
umn density range over which the cutoff is measured. This 
column density does in general not correspond to the mean 
gas density and hence logfeo will in general depend on 7. We 
will show below that this dependence can be removed by 
renormalising the equation of state. 

In Fig. M we plot the temperature at mean density pre- 
dicted from the power-law model (equation H) as a function 
of the true To. Data points are determined by fitting power- 
laws to 500 sets of 300 random absorption lines. Error bars 
enclose 68 per cent confidence intervals around the medians. 
Solid circles are used for data from simulated models, open 
squares are used for models created by imposing an equa- 
tion of state on model Ob. These conventions will be used 
throughout the paper. Fig. M is a similar plot for the slope of 
the equation of state (equation pi) . The predicted and true 
parameters of the equation of state are highly correlated. 
The slight offset between the predicted and true quantities 
simply reflects the fact that the optical depth weighted den- 
sity and temperature are not exactly the same as the true 
density and temperature of the absorbing gas (i.e. the T 
and p appearing in equation m are not exactly the same as 
those appearing in equations Hand H). The main conclusion 
to draw from these plots is that the power-law model works 
and that we can therefore use these equations to gain in- 
sight in the relationship between the equation of state and 
the cutoff in the 6(7V)-distribution. 

The objective is to establish the relations between the 
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Figure 7. The temperature at mean density predicted from the 
power-law model (equation tj) as a function of the true tempera- 
ture at redshift 2 = 3. Error bars enclose 68% confidence intervals 
around the medians, as determined from 500 sets of 300 random 
lines. Open squares indicate that the points resulted from impos- 
ing an equation of state with the given To on model Ob. All points 
have similar values of 7. The dashed line is the least-squares fit. 




Figure 8. The index 7 predicted from the power-law model 
(equation ra) as a function of the true index of the power-law 
equation of state at redshift 2 = 3. Error bars enclose 68% confi- 
dence intervals around the medians, as determined from 500 sets 
of 300 random lines. Open squares indicate that the points re- 
sulted from imposing an equation of state with the given 7 on 
model Ob, while leaving To unchanged. The dashed line is the 
least-squares fit. 
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Figure 9. The intercept of the b(N)-cuto5 as a function of the 
temperature at mean density at redshift z = 3. The error bars en- 
close 68% confidence intervals around the medians, as determined 
from 500 sets of 300 random lines. Open squares indicate that the 
points resulted from imposing an equation of state with the given 
To on model Ob. All points have similar values of 7. The dashed 
line is the maximum likelihood fit, as determined from the full 
probability distributions. The tight correlation and the relatively 
small errors make the intercept of the cutoff a sensitive measure 
of the amplitude of the equation of state. 



cutoff parameters and the equation of state using simula- 
tions. These relations can then be used to measure the equa- 
tion of state of the IGM using the observed cutoff in the 
6(iV)-distribution. The amplitudes of the power-law fits to 
the cutoff and the equation of state are plotted against each 
other in Fig. 0. The relation between log 60 and log To is lin- 
ear, implying that the coefficients appearing in equation H 
do not vary strongly with cosmology. The error bars, which 
indicate the dispersion in the cutoff of sets of 300 lines (typi- 
cal for z — 3), are small compared to the differences between 
the models. This means that measuring the cutoff in a single 
QSO spectrum can provide significant constraints on theo- 
retical models. 

The slope of the cutoff, V — 1, is plotted against 7 in 
Fig. hd. The relation between the two is linear, but F in- 
creases only slowly with 7. The dispersion in the slope of 
the cutoff for a fixed equation of state is comparable to the 
difference between the modelsfi. The weak dependence of F 
on 7 and the large spread in the measured V will make it 
difficult to put tight constraints on the slope of the equation 
of state. 



3.1 



Correlations 



The intercept of the cutoff in the fc(iV)-distribution is a mea- 
sure of the temperature at the characteristic density of the 
absorbers corresponding to the lines used to fit the cutoff. 
In general this is not the mean density and consequently 
the translation from the intercept bo, to the temperature at 
mean density To, depends on the slope 7. This is illustrated 
in the left panel of Fig. hll where the measured values of bo 



The simulated models all have similar values of 7 because they 
all have the same (/V-background. 



Figure 10. The slope of the 6(iV)-cutoff as a function of the in- 
dex of the power-law equation of state for redshift z = 3. The 
error bars enclose 68% confidence intervals around the medians, 
as determined from 500 sets of 300 random lines. Open squares in- 
dicate that the points resulted from imposing an equation of state 
with the given 7 on model Ob, while leaving To unchanged. The 
dashed line is the maximum likelihood fit, as determined from the 
full probability distributions. The slope of the cutoff is a measure 
of the slope of the equation of state, The weak dependence of T 
on 7 and the large variance in the measured T will make it hard 
to constrain 7 accurately. 



for models that have identical values of To, but a range of 7- 
values are compared. As predicted by the power-law model 
(equation 0), logfoo increases linearly with 7. 

In principle, the index 7 can be measured using the 
slope of the cutoff, to which it is proportional (equation pt). 
However, F increases only slowly with 7 and the statistical 
variance in V is large, making it hard to put tight constraints 
on 7. It appears therefore that even though log 60 is very 
sensitive to log To, and can be measured very precisely, the 
uncertainty in To will be relatively large due to the weak 
constraints on 7. It is important to realise that any statistic 
that is sensitive to the temperature, depends in general on 
both To and 7. It is for example impossible to determine To 
by fitting the ^-parameter distribution at N 7^ N(p). Since 
this statistic is sensitive to the temperature at the density 
corresponding to a column density N, any equation of state 
(To, 7) that has the correct temperature at this density will 
fit the data. 

Although it is conventional to normalise the equation 
of state to the temperature at the mean density, it can in 
principle be normalised at any density, ps = p(l + S) say, 

IogT = tagT a + (7-l)Iog(p/p«). (7) 

Similarly, equations H and ]4 can be generalised to 

log(p/p s ) = B\og(N/N s ), (8) 

log b = log 60 + (r - 1) log{N/N 6 ), (9) 

where N$ = N(p — p$) and consequently the coefficient A 
vanishes. In practice, the cutoff is fitted over a given column 
density interval and log Ng is set equal to the mean log N of 
lines in this interval. The measured cutoff is then converted 
into an equation of state normalised at the corresponding 
density, ps = p(log Ns = (log N)). Equation H then becomes 
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Figure 11. The intercept of the fe(A r )-cutoff as a function of the index 7. The models in the left panel all have the same To, while the models 
in the right panel have identical Tg, where 8 = 2. The intercept of the cutoff is a measure of the temperature of gas at the characteristic 
density corresponding to the absorption lines used to fit the cutoff (which is <5 ss 2 for the interval 10 12,5 cm -2 < N < 10 145 cm -2 at 
z = 3). Hence the intercept depends on the slope of the equation of state, unless it is normalised to the temperature at this density. 
Plotted data points are Ob (filled circle) and different equations of state imposed on Ob (open squares), all for 300 random lines at 
z = 3. Error bars enclose 68% confidence intervals around the median values. The dashed lines are the maximum likelihood fits. 



log 60 = (logTs — C)/D and the intercept of the cutoff de- 
pends only on the amplitude of the equation of state. 

At redshift z = 3, using the column density interval 
10 125 cm" 2 < N < 10 145 cm" 2 , N s = 10 136 cm" 2 and S « 
2. In the right panel of Fig. Ill] the intercept of the cutoff is 
plotted as a function of 7 for a set of models that all have the 
same temperature at this density. As expected, the intercept 
is insensitive to the slope of the equation of state. 

In summary, the intercept of the b(N)-cutoS is a mea- 
sure of the temperature of the gas responsible for the ab- 
sorption lines that are used to determine the cutoff. If we 
normalise the equation of state to the temperature at the 
characteristic density of the gas, then the intercept of the 
cutoff depends only on the amplitude of the equation of 
state. The slope of the cutoff is always determined by the 
slope of the equation of state. Hence the cutoff in the b(N)- 
distribution can be used to determine both Ts, where 5 is 
the density contrast corresponding to the mean log N of the 
lines used in the fit, and 7. The temperature at mean den- 
sity, To, depends on both Tg and 7 and therefore on both 
the intercept and the slope of the cutoff. 



4 MEASURING THE CUTOFF 

The main problem in measuring the cutoff in the b(N)- 
distribution is the fact that it is contaminated by spurious 
narrow lines. Line blending and blanketing, noise and the 
presence of unidentified metal lines all give rise to absorp- 
tion lines that are narrower than the lower limit to the line 
width set by the thermal state of the gas. We have there- 
fore developed an iterative procedure for fitting the cutoff 
that is insensitive to the presence of a small number of nar- 
row lines. In order to minimize the effects of outliers, robust 
least absolute deviation fits are used. The first step is to fit 
a power-law to the entire set of lines. Then the lines that 



have fa-parameters more than one mean absolute deviation 
above the fit are removed and a power-law is fitted to the 
remaining lines. These last two steps are repeated until con- 
vergence is achieved. Finally, the lines more than one mean 
absolute deviation below the fit are also taken out and the 
fit to the remaining lines is the measured cutoff. The algo- 
rithm works very well if there are not too many unphysically 
narrow lines. At the lowest column densities (^ 10 12 ' 5 cm -2 ) 
however, blends dominate and the cutoff is washed out. 

Fortunately, there are ways to take out most of these un- 
physically narrow lines. We have already shown (Fig. 0) that 
removing lines with large relative errors in the Voigt profile 
parameters significantly sharpens the cutoff. We choose to 
consider only those lines with relative errors smaller than 50 
per cent. A smaller maximum allowed error would result in 
the removal of many of the regular, isolated lines. 

Another cut in the set of absorption lines can be made 
on the basis of theoretical arguments. Assuming that ab- 
sorption lines arise from peaks in the optical depth r, and 
assuming that In r is a Gaussian random vari able (as is the 
case for linear fluctuations) , Hui & Rutledge ( 1997 ) derive a 
single parameter analytical expression for the fa-distribution: 



dN fa 4 
■« 0C F eXP 



b 4 



(10) 



where b a is determined by the average amplitude of the 
fluctuations and by the effective smoothing scale. Fig. |l2| 
shows the fa-parameter distribution for the lines plotted in 
Fig. and the best-fitting Hui-Rutledge function (dashed 
line). The b- value for which the Hui-Rutledge fit vanishes 
(dotted line) corresponds to the dashed line in Fig. H. The 
fa-distribution has a tail of narrow lines which is not present 
in the theoretical Hui-Rutledge function. These lines are in- 
dicated by diamonds in Fig. H. Direct inspection shows that 
all these lines occur in blends. Two examples are the lines in- 
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Figure 12. The 6-parameter distribution of the lines plotted in 
Fig. 0. The dashed curve is the best-fitting Hui-Rutledge func- 
tion (equation lOh for this distribution. The vertical dotted line 
indicates the point where this function has the value 10 -4 , which 
corresponds to the horizontal dashed line in Fig. H. Virtually all 
of the lines that are narrower than this cutoff in the Hui-Rutledge 
function are blends. We remove these lines from the sample before 
fitting the 6(JV)-cutoff. 



dicated by arrows in Fig. g. The size of the low 6-tail depends 
on the number of blended lines and therefore on the signal 
to noise of the spectrum and on the fitting procedure. How- 
ever, we find that in general, virtually all of the lines with 
&-values smaller than the cutoff in the fitted Hui-Rutledge 
function are blends. We therefore remove these lines before 
fitting the cutoff in the 6(A r )-distribution. 

Fig. h3 illustrates the effect of the two cuts (relative 
errors and Hui-Rutledge function). The probability distri- 
butions for the parameters of the fitted cutoff, bo and T, are 
plotted for different cuts. The dotted lines are the distribu- 
tions resulting from fitting the cutoff for the complete set 
of lines. Removing the lines with large errors or those with 
fa-values smaller than the cutoff of the Hui-Rutledge fit re- 
sults in a smaller intercept log b(Vl It makes no difference 
which cut is applied. Applying a cut in error-space does not 
affect the slope of the cutoff, T — 1. However, taking out 
the lines below the Hui-Rutledge cutoff, removes the low-b, 
\ow-N tail without affecting the higher column density end, 
and therefore yields a smaller slope. 

A typical QSO spectrum at z ~ 3 has about 300 Lya 
absorption lines between its Lya and Ly/3 emission lines 
with column densities in the range 10 12 ' 5 cm -2 < N < 



10 1 



' cm . The number density of Lya lines decreases 



rapidly with decreasing redshift. At z ~ 2 there are typi- 
cally less than 100 lines. The fact that the number of lines 
in an observed Lya forest is finite introduces statistical vari- 
ance. We therefore use many (500) realizations to determine 
the full probability distributions of the parameters of the 
cutoff. 

Fig. u4 illustrates the effects of changing the number 

T The removal of very narrow lines reduces the scatter around 
the cutoff and therefore causes the algorithm to converge at a 
lower mean absolute deviation. More iterations are needed before 
convergence is obtained, yielding a lower final cutoff. 
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Figure 14. The intercept (upper panel) and slope (bottom panel) 
of the b(7V)-cutoff for different numbers of absorption lines used 
in the fit. The cutoff depends weakly on the number of lines, 
although the statistical variance does increase significantly when 
the number of lines falls below 200. The numbers refer to the 
number of lines after removal of the lines with errors greater than 
50%, but before applying the Hui-Rutledge cut. The error bars 
indicate 68% confidence intervals around the medians. The lower 
curve in the bottom panel has been slightly displaced. 



of absorption lines. The algorithm is surprisingly insensi- 
tive to the number of lines. The parameters of the cutoff 
vary only slightly for 60 lines or more. The variance does 
decrease as the number of lines increases, but remains al- 
most the same for more than 200 lines. This suggests that 
the method should work even with the small number of lines 
per spectrum at z ~ 2. It also opens up the possibility of 
splitting higher redshift spectra into redshift bins. The fact 
that the cutoff depends weakly on the number of lines is not 
a problem, since we can determine the relation between the 
cutoff and the equation of state for any number of lines, in 
particular for the number of lines in an observed spectrum. 
While the statistical variance in the intercept of the 
6(A r )-cutoff is small, the variance in the slope is comparable 
to the differences between models. Fortunately, we can do 
better than measuring the cutoff for the complete set of lines 
in an observed spectrum. The bootstrap method (drawing n 
random lines from the complete set of n lines, with replace- 
ment) can be used to generate a large number of synthetic 
data sets. These data sets can then be used to obtain approx- 
imations to the probability distributions for the parameters 
of the cutoff. Since bootstrap resampling replaces a random 
fraction of the original lines by duplicated original lines, a 
smaller fraction of the b-N space around the cutoff is filled 



© 0000 RAS, MNRAS 000, 000-000 



Measuring the equation of state of the IGM 11 



0.10 
0.08 

0.06 

0.04 

0.02 

0.00 

1 







Ab/b.AN/N < 0.5 //,\t ■"■ 

■_ b >b min \ ■ -., : 

no constraints //, V ; 

: j \ \ : 


Ab/b,AN/N < 0.5 

b>b nin 

r ■ '■ no constraints 



1.2 



1.3 
log b (km/s) 



1.4 



0.1 



0.2 

r-i 



0.3 



0.4 



Figure 13. Probability distributions for the intercept (left panel) and slope (right panel) of the cutoff in the 6(A r )-distribution for 300 
random absorption lines from model Ob-hot at z = 3. The solid line is for the algorithm used throughout this paper. The other curves 
illustrate the effect of including lines with b-parameters smaller than the cutoff of the Hui-Rutledgc function (dashed), relative errors 
greater than 50% (dot-dashed) and both (dotted). All cutoffs are for the column density range 10 cm < N < 10 5 cm . 



and the variance in the measured cutoff increases. Although 
the bootstrap distribution is generally broader than the true 
distribution, its median is a robust estimate of the true me- 
dian. When dealing with observed spectra, we will use the 
medians of the bootstrap distributions as our best estimates 
of the parameters of the 6(Af)-cutoff. In section M we will 
use Monte Carlo simulations to estimate the variance in the 
bootstrap medians. 



5 SYSTEMATIC EFFECTS 

In section ra we established the relation between the cutoff in 
the 6(A jr )-distribution and the equation of state of the low- 
density gas. In this section we will investigate whether other 
processes can change this relation. 



5.1 Cosmology 

Cosmology affects not only the equation of state, b ut also 
determines the evolution of structure. Theuns et al. (1999a) 
showed that the 6-parameter distribution depends on cos- 
mology. In particular, they showed that the effect of peculiar 
velocity gradients on the line widths can be very different 
in different CDM variants. We have recomputed simulated 
spectra for model S after imposing the equation of state of 
the significantly hotter model Ob. In Fig. [15l the probabil- 
ity distributions of the parameters of the 6(A r )-cutoff are 
compared for model Ob (solid lines) and this new model, 
S-hot (dashed). The distributions are almost indistinguish- 
able. Also plotted is model Ob-vel (dot-dashed), which was 
created by setting all peculiar velocities in model Ob to zero. 
Again, the probability distributions are almost unchanged. 
The line widths of many of the low column density lines 
are dominated by the Hubble flow. Fig. Mm illustrates the 
effect that changing the Hubble expansion has on the cutoff 
in the 6(A)-distribution. The obvious way to change the 



Hubble flow is to change the Hubble constant. However, the 
equation of state also depends on the value of the Hubble 
constant. In order to isolate the effect of the Hubble flow, 
we changed the value of the Hubble constant in the analysis 
stage, i.e. just before computing the spectra, keeping the 
equation of state fixed. Increasing the value of the Hubble 
parameter at z — 3 from a corresponding present day value 
of h — 0.65 to 0.8 has no effect on the cutoff. Lowering h to 
0.5 shifts the slope to slightly larger values, but leaves the 
intercept unchanged. 

We conclude that, unlike the 6-distribution, the cutoff 
in the b(A)-distribution is independent of the assumed CDM 
model for a fixed equation of state. 



5.2 Mean absorption 

The optical depth in neutral hydrogen is proportional to the 
quantity nf/i 3 /rjjj. Since it is still unclear what the domi- 
nant source of the metagalactic ionizing background is, the 
Hi ionization rate, Tjjj, is uncertain. Cha nging r-rj] has 



very little effect on the equation of state ( Hui & Gnedin 
1997), which means that the optical depth can be scaled 
to match the observed mean flux decrement in the analy- 
sis stage. However, the observations do show some scatter 
in the mean flux decrement. This scatter could be due to 
spatial variations of the ionizing background, or it could be 
caused by measurement errors. We therefore need to check 
whether errors in the assumed effective optical depth affect 
the relation between the cutoff and the equation of state. 

Changing the effective optical depth by rescaling the 
ionizing background alters the relation between column den- 
sity and gas density. This will shift the 6(A r )-distribution 
along the A-axis. Hence we expect the intercept of the cut- 
off to change and the slope to remain constant. In terms 
of our power-law model, increasing the photoionization rate 
(i.e. decreasing the effective optical depth) will increase the 
coefficient A in equation ti and thus increase the measured 
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Figure 15. Probability distributions for the intercept (left panel) and slope (right panel) of the b(7V)-cutoff for models with identical 
equations of state, but different systematics. S-hot (dashed line) is model S with the equation of state of model Ob imposed. Ob-vel 
(dot-dashed) is identical to model Ob, but with all peculiar velocities set to zero. Ob+scat is model Ob with twice as much scatter 
around the equation of state. All distributions are determined using 500 sets of 300 random absorption lines, at z = 3. 
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Figure 16. Probability distributions for the intercept (left panel) and slope (right panel) of the b(7V)-cutoff for models with identical 
equations of state, but different Hubble flow contributions to the line widths. For models Ob-h=50 and Ob-h=80 the spectra have been 
calculated using a Hubble parameter corresponding to a present day value of h = 0.5 and h = 0.8 respectively (Ob has h = 0.65). All 
distributions are determined using 500 sets of 300 random absorption lines, at z = 3. 



intercept bo for a given To (equation Bh, while leaving the 
slope r — 1 constant (equation H). Fig. hjj confirms these 
predictions. The dependence of &o on the effective optical 
depth turns out to be rather weak, even for models with 
a relatively steep cutoff. Realistic errors in the mean flux 
decrement (10-20% at z ~ 3) will give rise to errors in bo 
that are smaller than the statistical variance. 



5.3 Signal to noise 

The signal to noise ratio (S/N) per pixel in the simulated 
spectra is 50, comparable to the noise level in high-quality 
observations. However, spectra taken with for example the 



HIRES spectrograph on the KECK telescope have a S/N 
that varies across the spectrum. Fig. Ilq illustrates the effect 
of changing the signal to noise (S/N) in the spectrum. The 
statistical variance in the cutoff increases rapidly when the 
S/N falls below 25. While the intercept increases slightly for 
S/N smaller than 25, the slope appears to be independent 
of the signal to noise. We therefore conclude that variations 
in the S/N in observed spectra are unimportant, as long as 
the signal to noise is greater than about 20. 
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Figure 17. The intercept of the fe(A r )-cutoff decreases slightly 
with increasing absorption (upper panel) . The slope of the cutoff 
and the effective optical depth appear to be uncorrelated (bot- 
tom panel). For realistic errors in the assumed effective optical 
depth (10—20% at z ~ 3), the error in the intercept is small com- 
pared to the statistical variance. Plotted data points are model Ob 
(filled circle) and models with rescaled background fluxes (open 
squares), all for 300 random lines at z = 3. All models have the 
same equation of state. Error bars enclose 68% confidence inter- 
vals around the median values. 



5.4 Missing physics 

Since the cosmological parameters and reionization history 
have yet to be fully constrained, the models used in this 
paper may not be correct. This will however not affect the 
results of this paper, as long as the bulk of the low-density 
gas follows a power-law equation of state. We should there- 
fore ask if there are any physical processes, which have not 
been incorporated in the simulations, that could destroy the 
uniformity of the thermal state of the low-density gas. Feed- 
back from massive stars and fluctuations in the spectrum 
of the photoionizing flux are two examples. Although these 
processes would predominantly affect the gas in virialized 
halos, they could introduce some additional scatter in the 
equation of state of the gas responsible for the low column 
density Lya forest. 

We recomputed simulated spectra for model Ob, after 
doubling the scatter around the fitted temperature-density 
relation. The fe(A r )-cutoff of this new model, Ob+scat, is 
indistinguishable from the one of model Ob (solid and dotted 
curves in Fig. |l5|) . Note that although processes like feedback 
can heat the IGM locally, it is hard to think of any process, 
apart from adiabatic expansion, that could cool the low- 
density gas. Furthermore, whereas low-density gas that is 



Figure 18. The intercept (top panel) and slope (bottom panel) of 
the &(7V)-cutoff as a function of signal to noise ratio (S/N). Vari- 
ations in the S/N are unimportant as long as the S/N is greater 
than about 20. Plotted data points are all for 300 random lines 
taken from model Ob at z = 3. All models have the same equa- 
tion of state. Error bars enclose 68% confidence intervals around 
the median values. 



heated to T > To(p/p) 7_1 can only cool over a Hubble time, 
gas that is cooled below the equation of state is quickly 
reheated. Any additional scatter in the thermal state is thus 
unlikely to alter the cutoff in the gas temperature. 

Any comparison between observations and simulations 
is complicated by the fact that observed spectra cover a 
much larger redshift path than the simulation boxes and 
therefore include effects like redshift evolution and cosmic 
variance, which are not present in the simulations. Large- 
scale fluctuations, especially in the ionization fraction (i.e. 
the relation between density and column density) could po- 
tentially distort the relation between the mean equation of 
state and the b(N)-cutoff. 

Redshift evolution mainly affects the mean absorption 
and does so in a well defined manner, which can be modeled 
by comparing the observations to a combination of simu- 
lated spectra that have different effective optical depths. If 
the line of sight to the quasar passes an ionizing source, 
then lines from that region will be shifted to lower column 
densities and therefore have little effect on the measured 
b(N)-cutoff (since it increases with column density). Since 
the ionizing background originates from a collection of point 
sources, minima in the ionizing background would be shal- 
lower and more extended than maxima, except during reion- 
ization. In any case, if the sightline goes through a region 
where the mean neutral hydrogen density is substantially 
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enhanced and therefore affects the measured cutoff, this will 
become clear when the spectrum is analyzed in redshift bins. 
We conclude that even if local effects would produce a 
large scatter in the thermal state of the low-density gas, the 
relation between the mean equation of state and the b(N)- 
cutoff would remain unchanged. 



6 MONTE CARLO SIMULATIONS 

In section W we described how bootstrap resampling can be 
used to reduce the statistical variance in the measured b(N)- 
cutoff. Given a set of absorption lines from a single QSO 
spectrum, the bootstrap method is used to generate syn- 
thetic data sets, for which the cutoff is measured. The me- 
dians of the resulting probability distributions for the pa- 
rameters of the cutoff are then used as best estimates of the 
true medians. 

In order to see how well the equation of state can 
be measured, we performed Monte Carlo simulations. We 
drew 300 random lines, with column density in the range 
10 12 ' 5 cm~ 2 < N < I0 145 cm- 2 , from model Ob at z = 3 
and used the bootstrap method to generate probability dis- 
tributions for the parameters of the cutoff. The median log bo 
and F were then converted to measurements of log To and 7 
using the linear relations determined from the simulations. 
The statistical variance in the medians of the bootstrap dis- 
tributions was estimated from 100 Monte Carlo simulations. 
For 300 lines per spectrum the dispersion in the median is: 
ffstat(logTo) = 0.033 (K), cr sta t( 7 ) = 0.13. 

If multiple spectra are available, the statistical variance 
in the median can be reduced by summing the bootstrap 
probability distributions of the different spectra. Fig. [13 il- 
lustrates how the errors change when multiple spectra are 
used, each containing 100 (dashed curves) or 300 (solid 
curves) absorption lines in the column density range over 
which the cutoff is fitted. The bottom curve of each pair 
indicates the statistical \a error. The statistical dispersion 
asymptotes to the bin size used for determining the proba- 
bility distributions. 

Besides the statistical variance in the median, there is 
some scatter around the linear relations between the param- 
eters of the cutoff and the equation of state. For example, 
for 300 lines per spectrum at z — 3, the systematic errors 
(i.e. the dispersion of the solid circles around the dashed 
lines in Figs. M and mh in the parameters of the cutoff are 
<T sys (log6 ) = 0.009 (kms~_M and cr sys (r) = 0.015. The top 
curve of each pair in Fig. ni| indicates the total (statistical 
plus systematic) la error. 

For a single spectrum with 300 lines (typical for z — 3), 
the predicted total errors are: cr(logTo) = 0.040 (K), (7(7) = 
0.17. For three spectra the errors reduce to cr(logTo) = 
0.029 (K) and 17(7) = 0.13. Adding more spectra has very 
little effect since systematic errors dominate. The errors are 
larger when there are only 100 absorption lines per spec- 
trum. In this case the statistical variance is larger and the 
errors can be reduced significantly by adding more spec- 
tra. The fact that even for this small number of lines the 
constraints on the parameters of the equation of state are 
significant, suggests that the method will also work at z ~ 2. 
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Figure 19. The dispersion in the amplitude (top panel) and slope 
(bottom panel) of the equation of state as a function of the num- 
ber of spectra used. The dispersion is computed from 100 Monte 
Carlo simulations of spectra with 100 (dashed curves) or 300 (solid 
curves) absorption lines each, at z = 3. The lower curves indicate 
the statistical variance, while the upper curves give the total er- 
rors (statistical and systematic). Combining spectra has only lim- 
ited effect since systematic errors dominate. The statistical errors 
asymptote to the bin size used for determining the probability 
distributions. 



7 SUMMARY AND DISCUSSION 

Numerical simulations indicate that the smooth, photoion- 
ized intergalactic medium (IGM) responsible for the low col- 
umn density Lya forest follows a well defined temperature- 
density relation. For densities around the cosmic mean, 
shock-heating is negligible and the equation of state of the 
gas is well-described by a power-law T — To(/?//5) 7_ . The 
equation of state depends on cosmology, reionization history 
and the hard X-ray background. Although the absorption 
spectra can be fitted by a set of Voigt profiles, the lines 
do not in general correspond to discrete high-density gas 
clouds. Scatter plots of the distribution of line widths (b- 
parameters) as a function of column density (N) in observed 
QSO spectra clearly show a lower envelope, which increases 
with column density. 

The decomposition of spectra produced by a fluctuating 
IGM into discrete Voigt profiles is artificial. However, the 
column density of the absorption lines correlates strongly 
with the density of the gas responsible for the absorption. 
Although the 6-parameters are in general not correlated with 
the temperature of the gas, the line widths of the subset of 
lines that are close to the 6(iV)-cutoff do show a strong cor- 
relation with temperature. This implies that there exists a 
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lower limit to the line width, set by the the thermal state 
of the absorbing gas, which in turn depends on its density. 
Hence the cutoff seen in the &(JV)-distribution is a direct 
consequence of the existence of a temperature-density rela- 
tion for the low-density gas and can be used to measure the 
equation of state of the IGM. 

We developed an iterative procedure for fitting a power- 
law, 6 = bo(N/No) T ~ 1 , to the &(iV)-cutoff over a certain 
column density range (10 125 cm" 2 < N < 10 145 cm" 2 at 
Z = 3). The algorithm is insensitive to unphysically narrow 
lines, which occur in blends and as unidentified metal lines — 



The intercept of the prtwpr-law, lngftn, can he measured very 
precisely and is shown to be very sensitive to log To (Fig. |9J). 
The slope of the cutoff, T— 1, is proportional to 7— 1, but the 
dependence is weak and it is harder to measure (Fig. |lOJ ), 

The intercept of the b(N)-cutoff is a measure of the tem- 
perature of the gas responsible for the absorption lines that 
are used to determine the cutoff. This gas is typically slightly 
overdense and consequently the intercept depends on both 
To, the temperature at mean density, and 7, the slope of 
the equation of state. However, if we normalise the equation 
of state to the temperature at the characteristic density of 
the gas, T = Ts(p/ ' ps) 1 " 1 , where ps = p(l + 8), then the 
intercept of the cutoff depends only on the amplitude of the 
equation of state, T$. 

The relation between the cutoff and the equation of 
state is independent of the assumed cosmology (for a fixed 
equation of state). In particular, it remains unchanged when 
all peculiar velocities are set to zero and when the contribu- 
tion of the Hubble flow to the line widths is varied. Changing 
the effective optical depth (i.e. rescaling the ionizing back- 
ground) , alters the relation between column density and gas 
density and thus the relation between log &o and log To . How- 
ever, the dependence is weak and realistic errors in the mea- 
sured mean absorption do not lead to significant errors in 
the derived value of To. Variations in the signal to noise ra- 
tio are also unimportant, as long as the ratio is greater than 
about 20. 

The simulations used to determine the relation between 
the cutoff and the equation of state do not incorporate some 
potentially important physical processes, like feedback from 
e.g. star formation. This will however not change the results 
presented in this paper, as long as the bulk of the low-density 
gas follows a power-law equation of state. Since local effects 
like feedback would increase the temperature and since gas 
cooled to a temperature below that given by the equation 
of state is quickly reheated, any additional scatter in the 
thermal state of the gas is unlikely to affect the cutoff in the 
gas temperature. Doubling the scatter around the equation 
of state has no discernible effect on the 6(A r )-cutoff. 

The finite number of absorption lines per QSO spectrum 
introduces statistical variance in the measured cutoff. The 
statistical variance can be reduced by using the bootstrap 
method to generate probability distributions for the parame- 
ters of the cutoff and using the medians as the best estimates 
of the true parameters. If multiple spectra are available, the 
variance can be further reduced by adding the bootstrap 
distributions of the different spectra. We use Monte Carlo 
simulations to estimate the statistical variance in the me- 
dians of the bootstrap distributions. Besides the statistical 
variance, there is a systematic uncertainty from the scatter 



in the linear relations between the parameters of the cutoff 
and the equation of state. 

For a single spectrum at z = 3 we predict the fol- 
lowing total (statistical plus systematic) errors: <r(logTo) = 
0.040 (K) (AT0/T0 = 0.09), cr( 7 ) = 0.17. For three spectra 
the errors reduce to cr(logTo) = 0.029 (K) (AT /T = 0.07) 
and 17(7) = 0.13. Increasing the number of spectra beyond 
three has very little effect on the total uncertainty because 
systematic errors dominate. These errors should be com- 
pared to the ranges considered to be physic ally reasonable, 
10 3.o K < To < IP 45 K and 1.2 < 7 < 1.7 ( |Hui, Gnedin& 
Zhang 1997] ; |Hui fc Gnedin 199l\ ) 



The analysis presented in this paper is for redshift 
Z — 3. At smaller redshifts, the constraints will be less tight 
because of the smaller number of lines per spectrum. How- 
ever, we showed that even for one third of the number of 
lines typical at z = 3, the constraints on the equation of state 
are significant. Furthermore, the statistical variance can be 
reduced by using multiple spectra. At higher redshifts the 
errors will also be somewhat larger than at z = 3, mainly 
because the higher density of lines increases the number of 
blends and the errors in the continuum fit. 

When using the simulations to convert the observed 
6(iV)-cutoff into an equation of state, one has to be careful 
to treat the observed and simulated spectra in the same way. 
For example, the simulated and observed 6(iV)-distributions 
should have the same number of lines and the same contin- 
uum and Voigt profile fitting algorithms should be used for 
the simulated and observed spectra. 

Given the existence of many high quality quasar absorp- 
tion line spectra, it should be possible to greatly reduce the 
uncertainty in the equation of state of the low-density gas 
over the range z = 2 — 4. This will allow us to put significant 
constraints on the reionization history of the universe. 
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